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CHAPTER 1: LANDING GEAR AERODYNAMIC NOISE PREDICTION 


USING UNSTRUCTURED GRIDS 

Frederic J. Souliez*, Lyle N. Long 1 ^ Philip J. Morris 1 and Anupam Sharma* 
Department of Aerospace Engineering 
The Pennsylvania State University 
University Park, PA 16802 

Abstract 

Aerodynamic noise from a landing gear in a uniform flow is computed using the 
Ffowcs Williams-Hawkings (FW-H) equation. The time accurate flow data on the 
surface is obtained using a finite volume flow solver on an unstructured grid. The 
Ffowcs Williams-Hawkings equation is solved using surface integrals over the landing 
gear surface and over a permeable surface away from the landing gear. Two geometric 
configurations are tested in order to assess the impact of two lateral struts on the sound 
level and directivity in the far-field. Predictions from the Ffowcs Williams-Hawkings 
code are compared with direct calculations by the flow solver at several observer 
locations inside the computational domain. The permeable Ffowcs Williams-Hawkings 
surface predictions match those of the flow solver in the near-field. Far-field noise 
calculations coincide for both integration surfaces. The increase in drag observed 
between the two landing gear configurations is reflected in the sound pressure level and 
directivity mainly in the streamwise direction. 
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D Wheel diameter 

H(f) Heaviside function, H(f) = 0 for/< 0 and H(f) = 1 for /> 0 
Li refer equation 2 

M local Mach number vector of the source 

M |M| 

Mr M i r i 

fi unit normal vector to the surface 

P~ compressive stress tensor 

p pressure 

p __ freestream pressure 

p' acoustic pressure (p- p„) 

Prms root mean square pressure perturbation 

r distance from source to observer 

r unit normal vector from source to observer 

ret retarded (source) time 

SPL Sound Pressure Level (dB) 

Tij Lighthill stress tensor 

t observer time 

U _ freestream velocity 

U i refer equation 2 

m i th fluid velocity component 

v,- i th surface velocity component of integration surface f = 0 
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x observer location vector 

8{f ) Dirac delta function, unity for /= 0, zero otherwise 

p density 

p 0 free-stream density 

T retarded time 

□ 2 wave operator 1 j c 2 d 2 /dt 2 -V 2 

Introduction 

The Ffowcs Williams-Hawkings (FW-H) equation has recently been used with 
permeable surfaces in order to predict aerodynamic noise 1,2 . It is an inexpensive method 
to include the quadrupole source terms inside the FW-H surface without performing any 
volume integrations. This can significantly improve the accuracy of the noise predictions 
at locations where nonlinear interactions in the flow cannot be ignored. This is notably 
the case with highly turbulent flows such as high Reynolds number jets and wakes. It is 
also only slightly more expensive to use than a moving Kirchhoff surface (see Ozyoriik 
and Long 3 ), but without the limitations of Kirchhoff methods. 

The motivation for predicting the far-field noise generated by a 4-wheel landing 
gear stems from the increasing contribution of airframe noise to the overall sound level of 
an aircraft in its landing approach. Early studies in the 1970’s by Heller and Dobrzynski 4 
showed that high-lift devices such as slats and flaps, as well as deployed gears, generated 
noise levels 10 dB higher than those of an aircraft in its “clean” cruise configuration. 
Aerospatiale (now EADS Airbus) investigated the noise produced by several Airbus 
airplanes, which seems to indicate that noise from high-lift devices is likely to dominate 
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for medium size aircraft, while landing gear noise seems more of a problem for existing 
and future high capacity aircraft. The importance of investigating landing gear noise is 
reinforced by Airbus Industry plans to extend the Airbus family towards a high capacity 
aircraft. 

Heller and Dobrzynski carried out a series of tests with both scale models 4 and 
full scale models 5 , which underscored the lack of detailed geometric features with model- 
scale experiments and their effects on high frequency noise. These early experiments 
also showed that there is an increase in noise radiation from tandem axle configurations, 
which is the second test case in the present study. However, it was also found during the 
full scale experiment that struts, braces and other small features contribute significantly 
to the overall sound level. A more recent work by Dobrzynski et al 6 where the impact of 
various gear sizes and configurations is measured, illustrates the difficulty in using scale- 
model results for full-scale noise predictions. The actual simulation of the landing gear 
flow field is also of interest since it potentially affects the inflow of flaps located 
downstream. This was experimentally shown by Stoker et al 7 during a wind tunnel 
investigation of the airframe noise radiated by a model-scale Boeing 777, in which case a 
second high-frequency noise source from the flap system is only seen in the presence of 
the landing gear. This landing gear - flap interaction noise source was even shown to 
increase significantly by using a highly detailed gear geometry. 

As already performed in a previous study by the same authors 8 , the goal here is to 
combine the flexibility of unstructured grids with the FW-H equation. We use the 
Parallel Unstructured Maritime Aerodynamics (PUMA) code for generating the flow 
data. PUMA has been validated in several instances for simulating time-accurate flow 
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data 9,10 . The aim in the present case is to evaluate the impact on the noise directivity and 
intensity of two landing gear geometries (LDG1 and LDG2). It is expected to observe 
larger pressure fluctuations and a more complex three-dimensional flow in the case 
involving two additional struts (LDG2). 

The Computational Grids 

The grids used for the simulation of the flow over both landing gear 
configurations were generated using the commercial package Gridgen by Pointwise, Inc. 
Figure 1 and Figure 2 show an overall view of the meshes on the landing gear surfaces 
with and without lateral struts respectively. The first mesh consists of about 80,000 
surface triangles, for a total of about 880,000 tetrahedra in the volume mesh. The second 
mesh reused as much of the previous grid features as possible. With two additional 
struts, the number of triangles on the surface went up to 135,000, with about 1.2 million 
tetrahedral cells. Specific attention was given to the cell clustering between the front and 
rear wheels, in order to capture as much of the wake from the upstream wheel impinging 
on the downstream wheel. Flow separation from the fore wheel and wake impingement 
on the aft are expected to generate large unsteady pressure fluctuations and therefore 
noise. With the second geometry, great care was given to the mesh refinement between 
the two lateral struts, with the aft strut in the wake of the fore strut. The smallest 
geometric features were not overly simplified, since they have been shown to generate 
high frequency noise as explained in a later section describing the Ffowcs Williams - 
Hawkings equation. 
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As shown in Figure 3 for the second geometry (LDG2), a porous FW-H 
integration surface was used in addition to the flow data collected on the landing gear 
surfaces themselves. This will help determine the magnitude of the quadrupole source 
term for this low Mach number flow. Permeable FW-H surfaces were used for both 
geometries, with about 13,000 triangles in the first case, and 15,500 triangles in the case 
including two lateral struts. This coarsening mesh away from the solid surface is due to 
computer limitations. It may not be able to support the higher frequency pressure 
fluctuations. The advantage of these porous surfaces is that they can capture quadrupole- 
like terms without having to perform any volume integration. FW-H surfaces can be 
used in regions dominated by nonlinear effects (unlike the Kirchhoff formulations). 

The Gibbs-Poole-Stockmeyer algorithm 11 was used to speed-up the 
communication process between CPUs. As shown in Figure 4, this procedure divides the 
domain into slices that minimize the number of messages between each processor, so that 
each CPU exchanges data with at most two neighboring CPUs. The time step for the 
unsteady simulations is determined by the smallest cell characteristic length. At a CFL 
number of 0.95, this yields a time step of 0.86E-08 second for the first grid, and 1 .90E-08 
second for the second grid (due mostly to some improved CAD work in the original 
geometry file). The numerical conditions were dictated by the CFL3D run performed at 
Langley: the Reynolds number based on the wheel diameter is 1 .25 million, for a free 
stream Mach number of 0.2. The actual wheel diameter of the landing gear model is 9.4 
cm. 
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Flow Solver 


PUMA (Parallel Unstructured Maritime Aerodynamics) is the computer program 
that was used to run the unsteady calculations. It is a finite volume, Runge-Kutta time- 
marching code that solves the compressible Navier-Stokes equations and uses 
unstructured grids. It uses the Message Passing Interface (MPI) library for parallel 
implementation. Its scaling performance for the two configurations is illustrated in figure 
5. The flop performance is slightly higher for the second case for any given number of 
CPUs, since the ratio of computation over communication is greater for a larger grid. 
The facility used to perform the computation is the latest our two Cost effective 
Computing Arrays (COCOA and COCOA2) 12 . COCOA2 is a Beowulf cluster comprised 
of 20 nodes each having dual 800 MHz Pentium III and 1 GB RAM. The cluster has dual 
fast-Ethemet on each node and all the nodes are connected using two HP2524 switches 
with channel bonding for increased data communication. These machines run Redhat 
Linux (version 7.0) and the gcc compiler. 

For these simulations the inherent artificial dissipation provided by Roe’s flux 
integration scheme acts as a sub-grid scale turbulence model. A parallel investigation on 
separated flow around a cone 13 shows that the implementation of a Large Eddy 
Simulation (LES) method using a Smagorinsky sub-grid scale model 14 may improve 
PUMA’s accuracy to simulate both mean and turbulent quantities in the wake of a cone 
base flow. LES has already been used extensively to compute sound sources 15 ’ 16 , but 
some recent work related to two-time statistics of LES data, would indicate that LES 
fields are too coherent if the eddy viscosity model does not include any random 
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backscatter 17 . One way to circumvent this may be the use of a dynamic LES, which is 
more likely to yield enough backscattering to decorrelate the fluid motion at large scales. 
An example of the use of a dynamic subgrid scale model combined with a Ffowcs 
Williams - Hawkings solver is given by Morris et al 18 in an attempt to simulate the jet 
noise for circular nozzles. 


Simulation Results 

Simulations were carried out over two cycles based on the expected shedding 
frequency of the wheel diameter. Each simulation took about 90 days on 24 CPUs. It is 
worth mentioning that the existing amount of data for both gear configurations put a 
strain on the available capacity in terms of storage requirements: about 40 Gigabytes of 
data have been collected for the two calculations described in this study. Data were 
sampled only for the second cycle to minimize the effects of the starting conditions. 
Local time stepping was initially used to accelerate the convergence from free-stream 
conditions to a realistic state. This is achieved by assigning to each cell the maximum 
allowable time step for a given CFL number (pseudo time marching). Global time 
stepping is then turned on before unsteady data is sampled. In order to evaluate the total 
drag, the momentum deficit method is used by evaluating the velocity deficit in the wake 
of the landing gear. More details can be found in Rae and Pope 19 . Figure 6 shows the 
average velocity deficit right behind the second gear configuration. In the second gear 
case, there is a good qualitative agreement with experimental results published by Stoker 7 
in the high-fidelity landing gear configuration. Figure 7 shows the drag forces computed 
by integrating the pressure on the gear surface (pressure drag) and using Pope’s wake 
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deficit approach (labeled total drag). Results are shown for both gear configurations 
during approximately 7 milliseconds of simulated flow time, giving enough time for the 
fluid to cover three times the gear total length. The force coefficients (drag, lateral and 
vertical force) are the computed forces divided by the landing gear surface area and the 
dynamic pressure. The increase in overall drag due to the introduction of the lateral 
support struts is large since these components are not aerodynamically profiled and are 
comparable to flat plates facing the incoming fluid flow. Figure 8 illustrates the lateral 
forces stemming from the presence of these two stmts. 

One expects the far-field sound pressure level to reflect this unsteady loading in 
both its intensity and directivity. Figure 9, which is a display of the instantaneous 
distribution of pressure on the landing gear surface in its second configuration, shows the 
pressure on these support stmts, as well as on the wheels. Figure 10 shows a 3D 
representation of some vortex filaments shedding off various gear components, and 
highlight the impact of the upstream elements’ wake onto gear elements at downstream 
locations. The effect of these vortices is not completely captured by the FW-H surface 
which lies on the landing gear itself. However the permeable FW-H surface does account 
for all the effects induced by these filaments until they cross its boundaries. 

Far-Field Noise Prediction 

Only recently has the FW-H equation been used on a permeable surface, di 
Francescantonio 20 was able to show that simply integrating the surface source terms on a 
porous FW-H surface does account for the quadmpole sources enclosed within the 
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surface. The FW-H equation is written in the standard differential form including all 
quadrupole, dipole and monopole source terms as 
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p'M 
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The subscript n indicates the projection of a vector quantity in the surface normal 
direction. Using the solution to the above equation given in Brentner and Farassat 21 and 
neglecting the quadrupole terms, the pressure fluctuation at a - given observer location x 
and time t is (equation 3 below) 
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FW-H code validation 

In the absence (at the moment) of experimental acoustic data to compare to, an 
already-proven method was implemented: the use of the CFD results to validate the FW- 
H sound predictions. As was done by the authors in a previous test case 8 , the pressure 
fluctuations computed by the flow solver PUMA at an observer in the near field were 
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compared with the predictions given by the FW-H post-processing utility. Although the 
near-field pressure fluctuations are large, and likely contain a great amount of 
hydrodynamic oscillations, the derivation of the FW-H equation is such that all pressure 
perturbations (acoustic and hydrodynamic) should be recovered. Examples in the near 
field are given by Farassat and Brentner 22 in the case of high-speed impulsive noise at 
rotor blade tip Mach number close to 0.9. It is assumed that at a Mach number of 0.2, the 
quadrupole terms do not contribute significantly to the far-field noise. The solution p Q 
to the quadrupole term of the FW-H equation is: 


4k p' Q (x,t) 


dx > dx j - /> o 


t J 

f f ^cdQdr 

j j r 


( 4 ) 


The volume integration, if performed, must be carried out over a large volume 
and represents a large computational task. The far field approximation of equation 4 
reduces to: 


= f f ^dQdr 

c to 2 1 f u r 


( 5 ) 


However, there is in the present case an interest in capturing quadrupole effects in 
the near field, so that an exact result to the FW-H is needed instead of the far field 
approximation. Farassat and Brentner 23 decomposed the quadrupole noise term into three 
components varying with 1/r, 1/r 2 and 1/r 3 respectively: 


- 11 - 



( 6 ) 



There is a possibility that the second and third terms may contribute in a 
significant way to the near field pressure variations. This implies that in order to validate 
the FW-H predictions against the CFD results one may have to account for some of these 
nonlinear effects in addition to loading noise in the near field since the observer is in a 
highly perturbed propagating medium. In the current derivation of the FW-H equation, 
the quadrupole term is not computed (to reduce computing time and to limit storage 
requirements). However the porous FW-H surface shown in a previous figure has the 
ability to recover all nonlinear effects occurring within its own boundaries. Figure 1 1 is 
an illustration of the instantaneous pressure distribution on the permeable FW-H 
integration surface. 

Observers were placed just above the landing gear main leg (x = 2.68 cm, y = 0 
cm and z = 17 cm for observer whose pressure is depicted in Figure 12), where the 
porous FW-H mesh is more refined, so that the FW-H predictions can take place using 
both the solid and the porous FW-H surfaces. An example of comparisons with the 
PUMA results at one of these near-field observer locations is shown in Figure 12. This 
shows good agreement between the porous FW-H surface predictions and the solver 
computation, whereas the solid FW-H surface misses by more than 50% some of the 
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pressure fluctuations. This demonstrates the ability of the second FW-H surface to 
predict the entire pressure oscillations, either of acoustic or hydrodynamic nature. It 
tends to suggest that quadrupole effects may represent a significant contribution to the 
overall near-field sound level even at moderate Mach numbers. Unless a volume 
integration is performed over the entire CFD domain, the entire pressure perturbation 
cannot be exactly reproduced where nonlinear effects are important and where vortices 
flow across the permeable FW-H surface. 

As expected, the agreement between the predictions from the two FW-H surfaces 
improves in the far-field. Figure 13 shows the pressure time history at 40 radii from the 
landing gear at a 50 degree angle with respect to the downstream axis. The field 
produced by the coarser FW-H surface off the landing gear does not reflect the same 
high-frequency fluctuations given by the predictions coming from data collected on the 
gear itself. In view of experimental results described earlier, it was decided to use the 
solid FW-H surface to investigate the far-field noise directivity, where high-frequency 
signals are thought to be significant. Figure 14 below illustrates the decrease of the RMS 
pressure signal as one moves away from the landing gear along the downstream axis. 
Calculations were made at 20 observers from 25 wheel diameters down to 50 wheel 
diameters in the wake of the landing gear. Both FW-H surface data were used and 
compared with a trend line assuming a signal decaying with 1/r. As observed previously, 
the agreement between the two surface predictions improves with increasing distance 
from the landing gear. As one moves further away from the landing gear, it is seen by the 
observer as an acoustic compact source, and the signal intensity should decrease with the 
inverse of the distance from the source. 
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Sound directivity patterns 

The sound directivity in the medium- and far-field requires the use of a parallel 
version of the FW-H post-processing program 8,24 . For each radii away from the landing 
gear, 72 observer locations are defined, so that a resolution of 5 degree angle is obtained. 
A total of 648 observer points were defined, and are illustrated in Figure 15, which shows 
the relative scale with respect to the landing gear. For both gear configurations, all three 
orientation planes were studied and the results are reported in Figures 16a to 16f, Figures 
17a to 17f and Figures 18a to 18f for the first and second gear configurations from a 
streamwise, spanwise and vertical perpective respectively. Polar directivity plots at 
radial locations of 10, 15 and 20 radii from the gear are plotted separately from the 
locations further away (30 to 50 radii from the gear) for scaling issues. Sound Pressure 
Level (SPL) contours with a reference pressure of 6x1 O' 5 Pa are also presented for radial 
locations varying from 25 to 50 radii from the landing gear. The scale for equivalent 
configurations is unchanged in order to allow for qualitative comparisons with respect to 
both directivity and intensity of the sound pressure signal. In Figures 16 to 18, the 
landing gear is not to scale, and is meant to illustrate which orientation axis is shown. 

The drag augmentation is reflected in the sound directivity patterns of both 
configurations. The intensity of the RMS pressure is greatly increased along the 
streamwise direction for the second gear case (LDG2). This is due to the two support 
struts on the aerodynamic profile of the landing gear. Regarding the lateral noise, the two 
struts seem to interfere with the build-up of sound, so that the signal in the spanwise 
direction is less than that in the clean configuration, where varying lateral forces on the 
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gear leg create pressure levels in the far-field comparable to those along the streamwise 
direction. In both cases, the near-field pressure perturbations are dominated by the 
fluctuations in drag. Little noise is generated in the vertical direction since most of the 
lift and drag variations are generated on the main gear leg. The overall pressure field 
looks much more disturbed in the second gear configuration, illustrating the complex 
three-dimensionality of the noise-generating flow pattern. 

Conclusion 

The flow field around two landing gear configurations of increasing complexity 
has been assessed. The wake deficit observed behind the landing gear is very similar to 
that experimentally measured on comparable configurations. A parallel version of the 
Ffowcs Williams — Hawkings equation has been implemented using inexpensive Beowulf 
clusters to extract near- and far-field sound information. Both solid and permeable FW-H 
integration surfaces have been used. Excellent agreement has been obtained in the near 
field between the porous FW-H surface predictions and the CFD solver results where 
hydrodynamic fluctuations are expected to dominate and are of greater magnitudes than 
those typical of acoustic signals. More work is needed in order to show that the 
discrepancy observed between the solid and porous FW-H surfaces in the near-field is 
linked to short-range quadrupole-like effects even though the problem that was dealt with 
in the present case is a relatively low Mach number flow for this kind of effects to be 
significant. 

The comparison of acoustic predictions produced by the two FW-H surfaces 
improves as the observer location is moved further away in the far field. The increase in 
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drag stemming from the lateral struts is reflected in the noise level and directivity. There 
is a significant increase in sound intensity in the streamwise direction, whereas the 
disturbance caused by these gear elements seems to interfere with the vortex shedding off 
the gear leg and the resulting lateral sound radiation. 
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Figure 3 Partial view of the porous FW-H surface around LDG2 gear configuration 








Figure 6 Average velocity deficit in the wake of the LDG2 configuration 



Drag Coefficient 



Figure 7 


Time-history of drag force coefficients for landing gear configurations 


LDG1 and LDG2 using pressure and wake methods 




Figure 8 Time-history of lateral (Fy) and vertical (Fz) pressure force coefficients 


for configurations LDG1 and LDG2 
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Figure 14 RMS pressure signal predictions from integration surfaces 1 and 2 and 


trend lines 
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Figure 15 Observer locations for sound directivity calculations 
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Figure 17 (top) Medium-field RMS pressure, (middle) far-field RMS pressure and 
(bottom) SPL contour from sideways perspective for LDG1 (left) and LDG2 (right) 





AERODYNAMIC NOISE PREDICTION 
USING PARALLEL METHODS ON 
UNSTRUCTURED GRIDS 

Lyle N. Long* Frederic Souliez* and Anupam Sharma 1 ’ 

Department of Aerospace Engineering 
The Pennsylvania State University , PA-16802 

Aerodynamic noise from a cone in a uniform flow is computed using the Ffowcs 
Williams-Hawkings (FW-H) equation. The time accurate flow data is obtained using a 
finite volume flow solver on an unstructured grid. The FW-H equation is solved for surface 
integrals over a permeable surface away from the cone. Predictions from the FW-H code 
are compared with direct calculations by the flow solver at a few observer locations inside 
the computational domain. A very good qualitative match is obtained. Sound directivity 
patterns in the azimuthal and in the longitudinal directions are presented. The FW-H 
code is also validated against a model problem of a monopole in a uniform mean flow. 


Nomenclature 

C p coefficient of pressure 

C 8 sub-grid scale constant in Smagorinsky model 

c sound speed in quiescent medium 

d base diameter of the cone 

f 8 vortex shedding frequency (Hz) 

H(f) Heaviside function, H(f) = 0 for / < 0 and 

H(f) = 1 for / > 0 
Li refer Eq. 2 

Lm LiMi 

L r Lifi 

L r L r fx 

M local Mach number vector of the source 

M |M| 

M n MiTli 

Mo Uq/c 

M r Mifi 

Mr Mifi 

n unit normal vector to the surface, n* 

Pa compressive stress tensor 

with p 0 Sij subtracted 
p pressure 

Po freestream pressure 

p ' acoustic pressure, p — p Q 

p’rms root mean squared pressure perturbation 

ret retarted time 

Tij Lighthill stress tensor 

t observer time 

0 angular location of the observer 

U averaged streamwise velocity 

Uq^Uoo freestream velocity 
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Uj refer Eq. 2 

Un UiUi 

Un UiTii 

Un UiTii 

Ui components of local fluid velocity 

u 'avg averaged streamwise perturbation velocity 

U n UiTii 

v n local normal velocity of the source surface 

6(f) Dirac delta function 

6(f) = 1 for / = 0, otherwise 6(f) — 0 

6ij Kronecker delta function, 

6ij = 1 for i = j, otherwise 6{j = 0 

p density of the fluid 

po freestream density of the fluid 

p f density perturbation, p - p 0 

fl vorticity (s -1 ) 

t o angular frequency of the monopole source 

□ 2 wave operator, d 2 ~ - V 2 ) 

Introduction 

R ECENTLY, the Ffowcs Williams-Hawkings 
(FW-H) equation has been used with permeable 
surfaces for predicting aerodynamic noise. The 
application of FW-H in this manner effectively allows 
for the inclusion of the quadrupole source terms inside 
the surface without performing volume integrations. 
This has significantly improved the accuracy of noise 
prediction for cases where the contribution from 
nonlinear interactions in the flow cannot be ignored. 
This is typical of highly turbulent flows, for example, 
high Reynolds number jets and wakes. 

The FW-H equation requires time accurate data on, 
and in the volume inside the permeable surface. This 
data is usually obtained by solving the Euler/Navier- 
Stokes equations accurately in time. Since the FW-H 
equation uses data from within the FW-H surface, the 
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outer grid can be made coarse without much loss of 
accuracy. Unstructured grids provide great flexibility 
in distributing the grid in the domain, and hence can 
be used to cluster the cells inside the FW-H surface. 
This feature can be exploited to significantly increase 
the computation speed while keeping almost the same 
accuracy in predicting aerodynamic noise. This will 
also permit the modeling of complex geometries such 
as helicopter fuselages, landing gear, and flaps. 

The goal here is to test the combination of unstruc- • 
tured grids with the FW-H equation in predicting the 
aerodynamic noise. The test case is chosen to be 
the flow over a cone. A cone has sharp edges which 
fixes the separation point. This makes the flow fairly 
Reynold’s number independent. 

We use the Parallel Unstructured Maritime Aero- 
dynamics (PUMA) 1 code for generating the time- 
accurate flow data. PUMA has been validated for 
time-accurate computations. 2 * ~ 4 The ultimate aim is 
to predict the airframe noise from complex geometries 
such as landing gear, slats, and flaps. This cone case 
may be considered as a benchmark problem. 

The Grid 

The grid used for the simulation of the flow over a 
cone of vertex angle 60° was generated using Gridgen. 
Figure 1 shows an overall view of the mesh consisting 
of approximately 280,000 tetrahedra. The clustering 
was done around the cone and in the wake region with 
increasing cell size towards the outer boundaries of the 
computational domain. The reason for using Gridgen 
comes from one interesting feature of this commer- 
cial software: arbitrary surfaces can be created around 
the cone (one within the CFD domain boundaries and 
the other being the CFD domain boundary) and are 
sources for the meshing algorithm. It is possible to 
export separately any of these closed surfaces in a sep- 
arate file, providing a means to extract flow data on the 
surface using a FW-H module that was added to the 
unstructured solver. The smallest cylinder was used as 
a porous FW-H surface. At the bounding faces of the 
CFD domain, Riemann boundary conditions were as- 
signed at each face center, hence minimizing reflections 
from the boundaries into the computational domain. 
The laxge cells in the far-field also help dissipate any 
reflections. A no-slip condition was used at the solid 
surface, even though the boundary layer was not re- 
solved due to computer limitations. 

By using a set of faces that are actually used by the 
flow solver during the computation, there is no addi- 
tional work required to extract the data needed for 
the far-field noise. This type of FW-H surface also re- 
flects the true mesh clustering present where the flow 
variables are locally being computed: there is no loss 
in accuracy due to the interpolation onto a surface 
whose refinement might not be that of the computa- 
tional grid. Since only the surface terms are evaluated 



Fig. 1 Overall view of the 280,000 cell mesh. 

during the acoustic prediction procedure, one does not 
have to take into account any phenomenon occurring 
outside the integration surface. The surface can also 
cross regions dominated by nonlinear effects. 

During a run, the faces (triangles in this case) would 
be identified and flagged on each CPU, so that face 
data would be output at a prescribed sampling rate 
(around 50 kHz in the present case): the sampling 
was done in such a'manner that one had at least 20 
data points per wavelength, the shortest wavelength 
being 10 times that of the simulated shedding fre- 
quency. To avoid any redundant data, faces shared 
between two adjacent CPUs had to be identified at 
the beginning of each run, so that the number of faces 
whose data are output is identical to the number of 
triangles on the actual FW-H surface. The grid par- 
titioning being done dynamically each time a run is 
initialized, the global cell indexing changes from run to 
run, making it necessary to run the above flagging pro- 
cedure any time the program is restarted. This makes 
the routine independent of the number of CPUs being 
used. Figure 2 illustrates the regions on the surface 
shared between 8 processors using the Gibbs-Poole- 
Stockmeyer reordering algorithm. 5 As expected, each 
region is a neighbor to at most two other partitions, 
minimizing the amount of inter-processor communica- 
tion. 

The time step needed for a time-accurate solution is 
determined by the smallest cell characteristic length. 
This is estimated to be one third of the cell volume 
divided by the maximum face area. For the grid 
described above, this yields a time step of 9.45E-08 
second at a CFL number of 0.9. The shedding fre- 
quency found during the experimental investigation of 
the flow is 36 Hz, for a Strouhal number equal to 0.171. 
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Fig. 2 Partitioning of the FW-H surface across 8 
processors. 

The Strouhal number was defined based on the cone di- 
ameter as St = f 8 d/Uoo. The numerical simulation is 
performed at Mach 0.2 at standard atmospheric pres- 
sure and temperature conditions, with an increased 
viscosity to match the experiment’s Reynolds number 
(50,000). Scaling the Strouhal number to the simu- 
lation’s Mach number yields a shedding frequency of 
230 Hz. The computation of a complete shedding cycle 
requires roughly 46,000 iterations. 

The Flow Solver - PUMA 

PUMA is a computer program, written in C, for 
the analysis of internal and external non-reacting 
compressible flows over arbitrary complex geometries. 
PUMA uses the Message Passing Interface (MPI) to 
rim the code in parallel. It can be run on arbitrary 
number of processors with very good scaling perfor- 
mance. Several papers 2,4 detail the benchmarking of 
the performance, and validation of PUMA. 

PUMA is based on finite volume methods and sup- 
ports mixed topology unstructured grids composed of 
tetrahedra, wedges, pyramids and hexahedra. The 
code may be run to preserve time accuracy for un- 
steady problems, or may be run using a pseudo- 
unsteady formulation to enhance the convergence to 
the steady state. Primitive flow quantities are com- 
puted at the cell centers. The code can be restarted 
from any point of time at which the solution is avail- 
able from previous computations. All flow variables 
are stored with double precision, but may be optionally 
stored as single precision to save memory and commu- 
nication time at the cost of reduced precision. 

Parallel Machines 

Computational Aeroacoustics (CAA) codes are usu- 
ally very computationally intensive. Even with very 


powerful machines, such jobs may require days, or even 
months to give results. Parallel computing using Be- 
owulf clusters offers an inexpensive way to handle such 
time-consuming simulations in reasonable amount of 
time. 

Three facilities offering parallel computational 
power at Penn State were used for the computations 
- COst effective Computing Array (COCOA), 2 CO- 
COA2 and LionX. 6 COCOA is a Beowulf cluster com- 
prising of 25 machines each having dual 400 MHz 
Pentium II processor. This facility was assembled by 
the authors and their colleagues in the Department of 
Aerospace Engineering at Penn State. The machines 
are connected via fast-Ethernet network which can 
support up to 100 Mbps bandwidth. A single Baynet- 
works 24-port fast-Ethernet switch with a backplane 
bandwidth of 2.5 Gbps is used for the networking. All 
the processors are dedicated to run parallel jobs. The 
operating system is Red Hat Linux. Message Pass- 
ing Interface (MPI) is used for parallel programming 
and the Gnu C compiler is used for compiling PUMA. 
Details regarding setting up and benchmarking of CO- 
COA may be obtained from Modi and Long 2 and 
COCOA’s website. 7 

COCOA was primarily set up to make parallel com- 
puting facility readily available to the CFD group of 
the Aerospace Engineering Department at Pennsylva- 
nia State University. The total cost of the cluster was 
just $80,000 in the year 1998, when it was set up. 
Since then this facility has been intensively used for 
various CFD simulations. COCOA2 is a newly assem- 
bled Beowulf cluster at Penn State. It has 21 nodes 
each having dual 800 MHz Pentium III processors and 
1 GB RAM each. The cluster has dual fast-Ethernet 
per node and all the nodes are connected using two 
HP2524 switches with channel bonding. 

Figure 4 plots the parallel speedup for COCOA and 
COCOA2 (1 Mflop = one million floating point opera- 
tions per second). Fairly good performance is obtained 
considering the small size of the problem. Figure 4 
shows the reduction in the flop rate per processor as 
the grid points are distributed over a larger number 
of processors. This trend is typical of Beowulf clus- 
ters as the ratio of computation over communication 
decreases. 

LionX is also a Beowulf cluster with 32 machines 
(each having dual 400 MHz Intel Xeon processors). 
These machines are connected via Myricom Myrinet 
with wire speed 1.28 Gbps. LionX also uses Linux with 
MPI for parallel programming. Performance com- 
parison and benchmarking results for LionX can be 
obtained from its website. 6 

CFD Results 

After initializing all variables to the freestream val- 
ues, local time stepping is used to accelerate the con- 
vergence towards a physically realistic flow. This is 
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Fig. 3 Parallel speed up for COCOA and COCOA2 



Fig. 4 

done by assigning to each cell the maximum allowable 
time step for a given CFL number based on each cell 
characteristic cell length. Global time stepping is then 
turned on for several cycles before data are sampled, 
to ensure that the data on the FW-H surface follow 
the equations of motion. Figure 5 illustrates the vor- 
ticity patterns in the wake of the cone, showing strong 
recirculation phenomena. The noise from this recir- 
culation is predicted by the FW-H module. Figure 6 
is the averaged streamline contour over one shedding 
period, illustrating the axisymmetric bubble that was 
observed during Calvert’s experimental study. 

In order to validate the solution, multiple compar- 
isons were made between the simulation and the exper- 
imental measurements. A basic Smagorinksy su -gn 
scale turbulence model 9 was added to the flow solver 
in order to improve the predictions, since a large-eddy 
simulation should yield better turbulent quantities. 



Fig. 6 Average streamlines over one shedding cy- 
cle. 

Figure 7 shows the averaged streamwise velocity 
profiles computed by the original flow solver and those 
computed by the same solver combined with an LES. 
In all three cases the magnitude of the reverse flow 
velocity is under predicted when compared with ex- 
perimental measurements. The predictions agree fair y 
well with Calvert’s data in terms of the length of the 
recirculation zone. Past the stagnation point, the re- 
sults including LES modeling follow the experimental 
curve more closely than those computed without any 
turbulence model. 

Figure 8 illustrates the variation of the pressure co- 
efficient C p along the wake centerline. In this case 
the LES having the largest sub-grid scale constant 
C. greatly over-corrects the pressure drop in the near 
wake of the cone. The LES using a Smagonnsky 
constant of 0.10 matches the measured pressure data 
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Fig. 7 Comparison of the averaged streamwise 
velocity with experiments for the flow solver with 
and without LES. 



Fig. 8 Comparison of the Cp coefficient with ex- 
periment for the flow solver with and without LES. 

very well until the stagnation point is reached. These 
results are consistent with those found in other re- 
lated investigations, using either the k-e turbulence 
model 10 or the k-e-v 2 model. 11 These simulations were 
compared against a set of experiments 12 at a lower 
Reynolds number (42,000). Madabhushi 13 also used 
an LES with as many as 850,000 mesh points, but 
completely over-predicted the length of the recircula- 
tion zone. 

Figure 9 shows that the averaged streamwise per- 
turbation velocity is not well predicted using any of 
the sub-grid scale constants. With the grid coarsening 
in the far wake, the fluctuating velocities are damped 
very rapidly as one goes away from the cone base. It 
is the flow solver without any turbulence model that 



Fig. 9 Comparison of averaged streamwise pertur- 
bation velocity with experiments with and without 
LES. 

yields unsteady velocity values that are closest to the 
experimental data. The solution without LES was se- 
lected to try to predict the fax-field noise. It also leads 
to the conclusion that a more advanced turbulence 
model (dynamic LES, Detached Eddy Simulation) is 
needed to simulate such separated flows, as found in 
Strelets. 14 

Far-Field Noise Prediction 

The two commonly used methods for far-field aero- 
dynamic noise predictions use the Kirchhoff equation 
or the Ffowcs Williams-Hawkings (FW-H) equation. 
While the governing equation in the ‘moving surface’ 
Kirchhoff formulation 15 is a convective wave equation, 
the FW-H equation is an exact rearrangement of the 
continuity and the momentum equations into the form 
of an inhomogeneous wave equation. Therein lies the 
strength of the FW-H equation over the Kirchhoff for- 
mulation. The FW-H equation gives accurate results 
even if the surface of integration lies in the nonlinear 
flow region. This is typically the case in jets and wakes 
when the nonlinear region extends to large distances 
downstream. 

In the Kirchhoff formulation the source terms are 
assumed to be distributed over a fictitious surface in 
the flow. The nonlinear effects (nonlinear wave prop- 
agation and steepening; variations in the local sound 
speed; and noise generated by shocks, vorticity, and 
turbulence in the flow field) happening within the 
Kirchhoff surface are captured by the surface integra- 
tion terms, but the Kirchhoff formulation requires the 
integration surface to be placed in a linear flow re- 
gion (i.e. far away from the body). This is difficult 
to achieve as most computational grids are generated 
with the concern of minimizing computations. Usually, 
a fine quality mesh is used near the body with increas- 
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ing cell size towards the outer boundaries. Therefore, 
the quality of the solution available in the linear flow 
region is generally bad. The FW-H equation, on the 
other hand, works fine even if the integration surface 
is in the nonlinear flow region. A detailed comparison 
of the Kirchhoff and FW-H formulations is provided 
in Brentner and Farassat. 16 

The solution of the full FW-H equation requires the 
evaluation of two surface integrals and one volume 
integral. The surface integrations correspond to the 
“thickness” noise (monopole) and the “loading” noise 
(dipole). The volume integration corresponds to the 
quadrupole term which accounts for the nonlinearity 
in the flow. 15,17 Evaluating the volume integral can 
be extremely computationally intensive and difficult 
to implement. Fortunately, the quadrupole term can 
be safely ignored for most subsonic flows as is the case 
in the present study. 

Only recently has the FW-H equation been used on 
a fictitious (i.e. not the same as the body) permeable 
integration surface 18 - exactly like the Kirchhoff ap- 
proach. di FYancescantonio 18 demonstrated that when 
the FW-H approach is applied on a Kirchhoff- type 
surface, the quadrupole sources enclosed within the 
surface cure accounted for by the surface sources. It 
should be noted that the “thickness” noise and the 
“loading” noise as obtained from solving FW-H equa- 
tion do not have any physical significance if the surface 
of integration is chosen to be permeable (fictitious). 
However, when the integration surface coincides with 
the body, these terms provide a physical insight into 
the source of sound generation. 

The FW-H equation is written in the standard dif- 
ferential form as 

°V(M) - asbiw <» 

- £:[Li6(f)} + j t [(poU n )6(f)] 

where Li and U n are defined as 

U n = Uihi :: Ui = (1 - —)vi + — 

Po Po 

Li = Pijrij + pUi{u n - v n ) (2) 

and Tij is the Lighthill stress tensor. The FW-H equa- 
tion can be solved using the formulation in Brentner 
and Farassat, 16 and the solution can be written in an 
integral form as 
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Fig. 10 Validation of the FW-H code against the 
analytical solution for a stationary monopole in a 
uniform mean flow. 


The quadrupole term is ignored in the present formu- 
lation. The integrations are performed on the FW-H 
surface at retarded time. Since the FW-H surface is 
fixed relative to the body (the cone) for this study, and 
the flow Mach number is constant, the following terms 
in the above integrals are zero : Un = M r = 0. The 
standard time binning technique discussed by Ozyoriik 
and Long 19 is used for obtaining pressure at the ob- 
server locations. 

The FW-H code and its Validation 

The FW-H code is written in Fortran 90. The 
code was tested for a model problem - a stationary 
monopole in a uniform mean flow. The FW-H surface 
is chosen to be a box made up of rectangular panels. 
The analytical solution to the model problem is eval- 
uated at the center of each panel to obtain the time 
history of the primitive variables on the FW-H sur- 
face. The prediction from the FW-H code (using the 
analytical data on the surface as input) is then com- 
pared with the analytical pressure perturbation at a 
point outside the surface. Figure 10 compares at an 
arbitrary point (300 m, 0, 0) the pressure perturba- 
tion predicted by the FW-H code and that obtained 
analytically for a stationary monopole source with an 
amplitude of 0.01 Pascals and a frequency of 2.267 Hz 
placed in a uniform mean flow of 0.3 Mach number. 
The analytical solution to this problem is : 


0(x,f) = 


e exp (twr.) 


47r [(x + Uo(t* - t)) 2 + y 2 + z 2 ] 1 / 2 

1 

i . Mo(x+Up(r» -t)) 
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[(i+Uo(r,-t)) 2 +i/ 2 +x 2 ] 1/2 


where r* is given by 
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Fig. 11 Comparison of the FW-H prediction us- 
ing unstructured surface grid against the analytical 
solution. 


The unstructured grid over the cone is created such 
that there is an unstructured cylindrical surface en- 
closed in the computational domain (Fig. 1). This 
surface is chosen to be the permeable FW-H surface. 
The elements of the surface are faces of the tetrahedra, 
and therefore, triangles. Since these triangles are cho- 
sen from the unstructured mesh, the area and normal 
varies from element to element. This, however, is not 
a problem because the FW-H equation only requires 
information on a closed surface; it does not depend on 
the structure of the elements constituting the surface. 
Clustering of the surface elements is desired to increase 
the resolution of the sources. The FW-H surface used 
for the present computation is the inner cylinder in 
Fig. 1. This grid was used with the model problem of 
stationary monopole in a uniform mean flow to test if 
the unstructured grid poses any problems. A perfect 
match is observed between the FW-H prediction and 
the analytical solution (Fig. 11). The comparison is 
made at an aribtrary point (300 m, 0, 0). This con- 
firms that an unstructured-mesh surface can be used 
as a FW-H surface without any loss of accuracy. Note 
that the first few seconds where the FW-H prediction 
does not match the analytical solution is the time it 
takes for the sound to reach the observer. This delay 
is more in Fig. 11 than in Fig. 10 because the unstruc- 
tured FW-H surface is very small and hence, farther 
away from the observer point than the structured sur- 
face used for Fig. 10. 

Results for the Cone 

PUMA is used to obtain time accurate data (the 
primitive flow variables) on the FW-H surface. One 
complete shedding cycle of the simulation is used for 
far-field noise prediction. Pressure at a few points out- 
side the FW-H surface (in the near field) is collected 
to compare with the predictions of the FW-H code. 
Four points distributed in the azimuthal direction near 
the base of the cone and very close to but outside the 
FW-H surface were chosen for comparison. The co- 


Point No. 

x (m) 

y (m) 

z (m) 

1 

0 

-0.055 

0 

2 

-0.025 

0.055 

0 

3 

0 

0 

0.05 

4 

0 

0 

-0.055 


Table 1 Coordinates of the observer locations for 
comparing FW-H predictions against PUMA. 

ordinates of the points are tabulated in Table 1. The 
cone has a base diameter of 0.02 m and a vertex angle 
of 60°. The center of the base of the cone is at the ori- 
gin and the vertex points upstream (positive x). The 
FW-H surface is a cylinder of radius 0.05 m and length 
0.175 m, centered at the origin. 

Figure 12 compares the pressure fluctuations at the 
four points listed in Table 1. Note that the PUMA 
pressure predictions have been shifted up by 20 Pas- 
cals. This is relatively a very small amount, about 
0.02% of the mean pressure. We believe that this 
under-prediction by PUMA may be due to the dis- 
sipation caused by inadequate clustering of grid cells. 
It may also be due to the small sample size, and we 
plan to do ensemble averaging. Note that this error 
is of the order of magnitude of pressure pertubations 
predicted by the FW-H code at any point inside the 
FW-H surface, which should actually be zero. How- 
ever, the prediction by the FW-H code agrees very well 
qualitatively with the PUMA solution. 

Sound Directivity 

The directivity of the noise from the cone was ob- 
tained by calculating the root mean squared (r.m.s.) 
pressure perturbation for one shedding cycle at differ- 
ent observer locations in azimuthal and longitudinal 
directions. Since the calculation for one observer loca- 
tion is completely independent of any other location, it 
is a perfect problem to run in parallel. Long and Brent- 
ner 20 suggested some self-scheduling parallel methods 
for multiple serial codes. However, no parallalization 
was done for the noise prediction results presented 
here. 

Figure 13 plots the directivity pattern in the az- 
imuthal direction on the plane x = —0.1 m, which is 
right behind the base of the cone. The pattern in Fig. 
13 is symmetric because of the symmetry of the cone 
about its axis. Since the FW-H equation cannot pre- 
dict the pressure fluctuation inside the FW-H surface, 
we can compute the noise only outside the FW-H sur- 
face. Therefore, the directivity patterns cure plotted in 
an annular region outside the FW-H surface. 

Figure 14 plots the directivity pattern in the longi- 
tudinal direction on the z = 0 plane. Since the noise 
is caused by both turbulence and fluctuating surface 
forces, the directivity shows several lobes. 

A conventional polar directivity pattern in the lon- 
gitudinal direction (z = 0 plane) is plotted in Fig. 15 
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Fig. 12 Comparison of pressure fluctuation, p — poo 
as predicted by PUMA and FW-H code at various 
locations listed in Table 1. 8 



Fig. 13 Directivity of the noise in the azimuthal 
direction behind the base of the cone ( x = —0.1). 



Fig. 14 Directivity of the noise in the longitudinal 
direction in the plane z = 0. 

for observers at 10 different radial locations (r = 0.15 
- 0.24 m). In Fig. 15, the cone is pointing to the right; 
the radial distance from the origin is equal to the r.m.s. 
pressure and the angle (theta) illustrates the location 
of the observer point in the domain. 

Conclusions 

Aerodynamic noise from a cone has been studied 
as a model problem to test the possibility of using 
unstructured grids for noise prediction from compli- 
cated bodies like landing gears, slats etc. A finite 
volume flow solver, PUMA has been used to obtain 
time-accurate flow data on a permeable FW-H sur- 
face. The FW-H code was validated against a model 
problem of a monopole in a uniform mean flow. The 
predictions from the FW-H code have been compared 
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Fig. 15 Polar plot of sound directivity in z — 0 
plane at a few radial locations. 


at four observer locations in the near field with direct 
calculations from PUMA. Noise predictions are made 
for a period of one shedding cycle. The comparison 
is fairly accurate with only a small D.C shift error. 
The directivity patterns of the noise from the cone 
are plotted in azimuthal and longitudinal directions. 
The sound directivity pattern has been shown to be 
fairly complicated due to the complex physics inside 
the FW-H surface. 
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Abstract 


The wake-dominated unsteady subsonic flow of Mach number 0.2 past a cone of vertex an- 
gle 60" is calculated numerically ustng high-order finite difference schemes on structured grids. 
The three-dimensional compressible Euler equations are solved to simulate an inviscid flow 
tha, exhibits large fluctuations of pressure and velocity due to the sheddtng of vortices behind 
fte cone. An axisymme.ric shuctured grid system is used, that is generated by rotating a two- 
dimensional grid plane around a centerline. Tire grrd singularity a, the centerline, where the 
Jacobian and some grid metrics approach infinity, is avoided by changing the form of the flux 
vectors in the Euler equations without any asymptotic assumption or simplification. Fourth- 
and sixth-order finite difference schemes are used for the evaluation of spatial derivatives, and 
a fourth-order Runge-Kutta scheme is used for marching the solution in time. The complex 
wake motions behind the cone are investigated by visualising the vorticity field. The mean flow 
pattern and periodic phenomena arc analyzed and compared with the experimental data. This 
demonshates the accuracy of the present approach to further analyses of wake-domma.ed 
flows past axisymmetric bluff bodies. 


I. Introduction 

Calvert [1] studied experimentally the low-speed flow patterns and the associated periodic 
phenomena with flow pas, a cone of various vertex angles. He noted that there was little pub- 
lished work on incompressible flow pas. an axisymmetric blunt-based body, except some ex- 
periments carried out between the 1930's and 1960's. These earlier studies went restricted to 
flow visualization a, low Reynolds numbers. He also discussed Are periodic phenomena associ- 
ated with such flows that are presumably related ,0 some kind of regular vortex shedding: 
however, the actual wake pattern is unknown above a Reynolds number of a few hundred. Cal- 
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vert's measurements of mean and fluctuating velocity, and mean pressure for various cone an- 
gles show that the wakes are all similar. This problem has received little attention in either ex- 
perimental or computational studies during the last three decades. Recently, Long et al [2] per- 
formed a simulation of a cone wake flow using low-order finite volume methods on unstruc- 
tured grids. The wake pattern and periodic phenomena in subsonic flow past a blunt-based 
body are of importance in the design of airframe components such as landing gears, especially 
for the reduction of aerodynamic noise: but, the existing knowledge of their flow characteristics 
is still far from complete. This paper seeks to provide some understanding of these complex 
flows through the simulation of the wake generated by an axisymmetric blunt-based cone. 

The present paper focuses on the application of high-order finite difference schemes and 
structured grids to the computation of a three-dimensional, unsteady, inviscid, compressible 
flow past a cone at subsonic speeds. The scope of the present study is limited to the large scales 
of the wake turbulence and the shedding vortices, rather than the fine scales of turbulence 
whose simulation would demand huge computational resources and time. Therefore, the com- 
pressible Euler equations are used in the present computation. Fourth- and sixth-order numeri- 
cal schemes, based on central finite differences, are used for the solution of the Euler equations. 
These methods have been developed for high accuracy and high resolution in the direct compu- 
tation of unsteady flows and acoustics [3, 4]. In addition, an artificial dissipation model and 
time-dependent boundary conditions are combined with the high-order schemes for a long- 
time stable solution [5-7]. 

The solutions are produced in a structured grid system generated by rotating a two- 
dimensional grid plane around the axis of symmetry. The axisymmetric structured grid has a 
grid singularity at the centerline, where the Jacobian and some grid metrics approach infinity. 
When a conventional form of the Euler equations is considered, this singularity makes it hard to 
calculate the flux variables at the centerline and the flux derivatives near the centerline when 
wide differencing stencils are used across the centerline in a generalized coordinate system. 
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This paper proposes a way to avoid the centerline singularity by changing the form of the flux 
vectors in the Euler equations without any asymptotic assumption or simplification. In this 
way, all the flux variables and derivatives may be evaluated without loss of accuracy. Tire sin- 
gularity treatment is a great help in the calculation of the flux derivatives near the centerline: 
however, the Euler equations are still indeterminate at the centerline. In the present approach, 
the solutions themselves at the centerline are interpolated from the neighboring values. 

The problem conditions for the present computation are a free-stream Mach number of 0.2 
and cone vertex angle of 60°. The present work provides a visualization of the wake flow field 
with its shedding vortices, and an analysis of the mean flow patterns and the periodic phenom- 
ena. The results are compared to Calvert's experimental data [1] for validation of the accuracy 
of the solution. A spiral pattern of the wake flow behind the cone is shown clearly through a 
visualization of the wake vorticity field. The initial location of the regular vortex shedding and 
its Strouhal number are simulated correctly. These successful comparisons indicate that the pre- 
sent methodology is capable of describing three-dimensional, unsteady, bluff-body wake flows. 

In the next section, the governing equations and the treatment of the centerline singularity 
are described. The high-order finite difference schemes are then introduced. The grid, and ini- 
tial and boundary conditions are then discussed. Both qualitative and quantitative results for 
the wake flow field and a comparison with experiment are then presented. 


II. Governing Equations and Centerline Singularity 

The governing equations are the three-dimensional compressible Euler equations. The flux 
vector form of the Euler equations transformed to the computational domain may be expressed 
in generalized coordinates as. 


3Q 3E 3F 3G A 
— + — + — + = 0 

dt d£ dv dC 
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where the flux vectors in generalized coordinates may be represented as 



E = y (£E + £F + £G), F = - fo,E + 77 y F + 77,0), 


G=y(^E + C y F + ^G). 


( 2 ) 


The conservative variables and the flux vectors in Cartesian coordinates are given by 
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where the total energy per unit mass is defined as e, =p/[(y-l)p]+(w 2 +v 2 + w 2 )/2 and 
y = c p jc v is the ratio of specific heats. y-\A in the present computation. In the Eq. (2), the 
transformation Jacobian, / and the grid metrics, are given by 
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In an axisymmetric structured grid, the Jacobian approaches infinity at the centerline and 
the denominator in Eq. (3) becomes zero. Therefore, all the grid metrics expressed in Eq. (4) are 
indeterminate at the centerline. However, one may easily show that the ratios of the grid met- 
rics to the Jacobian remain finite or zero at the centerline, even though the Jacobian and some 
grid metrics themselves have infinite values. This suggests a way to remove the centerline sin- 
gularity. The key is to choose a different way of expressing the Jacobian and the grid metrics. 
New variables that replace the conventional Jacobian and the grid metrics are defined by 
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( 5 ) 


( 6 ) 


The new variables with a superscript asterisk defined in Eqs. (5) and (6) have finite or zero val- 
ues in the entire computational domain, including at the centerline. Using the new variables, 
the expressions for the flux vectors given by Eq. (2) become 

Q = J'Q, E=£E+£F+£G, F = 77;e + 77;f + 77;G, G = £E + £F + £G. (7) 


In this manner, the centerline singularity is removed and all the flux variables at the centerline 
may be evaluated. Especially, this treatment of the singularity benefits the calculation of the flux 
derivatives in Eq. (1) near the centerline when using wide differencing stencils across the center- 
line. However, the Euler equations still cannot be solved at the centerline since J‘ is zero there 
and the vector of the conservative variables Q vanishes from Eq. (7). This means that Eq. (1) 
becomes indeterminate and it is impossible to integrate the solutions in time. However, the so- 
lutions (conservative variables) at the centerline may be interpolated from the neighboring val- 
ues. 


III. High-Order Finite Difference Schemes 

High-order finite difference schemes with high-resolution characteristics are used in the 
present computation on a structured grid. The main scheme is a pentadiagonal type of central 
compact finite difference scheme [3, 4). It is a generalization of the seven-point stencil Pade 
scheme used on the interior nodes. It may be expressed as 
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(8) 
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where f. is an objective function for the flux variables and f' is its spatial derivative on the i- 
th node. The grid spacing h is a constant independent of the index i in the computational do- 
main where all the grid points are equally spaced. Equation (8) may be solved by inverting a 
pentadiagonal matrix. The matrix must be completed at the boundaries. Therefore, non-central 
or one-sided formulations other than Eq. (8) are needed on the boundary and the near- 
boundary nodes to complete the matrix. These may be expressed as 


1 3 

/o + «WT+ Poafi = T Z a o* (/« ” fo ) for z==0 (boundary node), 
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(10) 


( 11 ) 


The coefficients in Eqs. (8)-(ll) are listed in Table 1. They are optimized, as described in Refs. 3 
and 4, to achieve maximum resolution characteristics with fourth-order accuracy (second-order 
in Eq. (9) for numerical stability). The optimized fourth-order compact finite difference schemes 
are used to evaluate the flux derivatives in the mean flow and radial directions. On the other 
hand, the flux derivatives in the azimuthal direction are calculated with a conventional central 
finite difference scheme for a simple handling of the periodic condition across the branch cut. 
The standard central finite difference scheme is given by 


45 (A, zlnhMii zLihlh* zM 

J> 6 Oh 


(12) 


Equation (12) has sixth-order accuracy, which is the highest order of truncation for the given 


6 



standard seven-point stencil. Combined with the high-order finite difference schemes in space, 
the classical fourth-order four-stage Runge-Kutta scheme is used for marching the solutions in 
time. 

As mentioned in Section II, the solutions (conservative variables) are obtained at the center- 
line not by solving the governing equations but by interpolation from the neighboring solu- 
tions. In the present work, a fourth-order interpolation is used at the centerline. It may be writ- 
ten as 

h ~ (13) 

where the negative indexes mean the values in the opposite direction across the centerline. Ac- 
tually, there are as many sets of the centerline solutions interpolated by Eq. (13) as half the 
number of grid points in the azimuthal direction. A unique value of the centerline solution is 
finally acquired by averaging these values. 

High-order schemes in space and time resolve a wider range of wavenumber or frequency 
than low-order methods. However, even the present schemes do not resolve the high 
wavenumber or frequency range effectively, an adaptive nonlinear artificial dissipation model 
[6] is also used to remove the unwanted numerical oscillations that may develop from the unre- 
solved range. The artificial dissipation model is implemented only at the last (fourth) stage of 
the Runge-Kutta scheme in order to minimize computational costs. In addition to the stringent 
requirements on the high-order and high-resolution numerical schemes, an accurate and robust 
calculation depends heavily on the suppression of any waves that may result from unwanted 
reflections on computational boundaries. The boundary conditions for a time-dependent prob- 
lem should be physically correct and numerically well posed. Generalized characteristic bound- 
ary conditions [7] are used as the time-dependent boundary conditions in the present computa- 
tion. Non-reflecting inflow/ outflow and the inviscid wall conditions are imposed at the 
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boundaries of the computational domain and the cone surface respectively. 

The accuracy of the high-order compact finite difference schemes, the adaptive nonlinear 
artificial dissipation model and the generalized characteristic boundary conditions has been 
validated through the previous publications of Refs. 3 to 7. They were applied to multi- 
dimensional steady/ unsteady and inviscid/ viscous computations in the various benchmark 
problems: linear wave convection in Ref. 3, acoustic radiation from axisymmetric baffle in Ref. 
4, noise generation from airfoils in Ref. 5, shock-sound interaction in a transonic nozzle in Ref. 6, 
and wake flow and acoustic radiation from a circular cylinder in Ref. 7. It was shown that they 
produced very accurate numerical solutions in comparison with analytic solutions and experi- 
mental data. They can be used effectively in the present computation. 

IV. Procedure and Results of Numerical Simulation 

In this section, the numerical simulation of an unsteady inviscid compressible flow past a 
cone of vertex angle 60° in a subsonic speed of Mach number 0.2 is described. The simulation 
uses high-order finite difference schemes and the modified form of the flux vectors in the Euler 
equations, in order to eliminate the centerline singularity. 

A. Grid System and Computation Procedure 

For the present problem, the computational domain consists of two blocks and rotating a 
two-dimensional grid plane around the centerline as shown in Figs. 1 and 2 produces an axi- 
symmetric grid system. The number of grid planes in the azimuthal direction is 48. The number 
of grid points is 101x50x48 in block 1 and 161x51x48 in block 2. This gives a total of 636,528 grid 
points. The grid points are clustered along the centerline and near the edge of the cone base. 
The minimum grid size at the center of the base is Ax/D = 0.007, Ay/D = 0.007 and Az/D = 
27tx0.007/48 where D is the base diameter. 
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The time step size is determined by the CFL (Courant-Friedrichs-Lewy) condition [8] which 
is given by. 


At = 


min 


ChJ' 


i^i+n+rMvc < W +^; 2 +vc +c +c y 


(14) 


with. 


U' = £« + £v + , F’ = + ^‘v + rj]w, 


w-m£u+c;v+c;w 


where the Courant number C = 1.0 is used in the present computation and c is the local speed of 
sound. In the evaluation of the time step, the operator 'min' in Eq. (14) does not include the cen- 
terline as the denominator approaches infinity there and the Euler equations are not solved on 
the centerline as remarked in the previous sections. 

A steady-state flow field is first acquired by solving the two-dimensional axisymmetric 
Euler equations. This is used as the initial condition for the present three-dimensional computa- 
tion. Figure 3(a) shows the pressure contours and Fig. 3(b) shows the Mach number contours. 
There is a recirculation bubble in cone base region. Calvert [1] deduced this characteristic form 
from his experimental measurements. The calculation, using the OpenMP [9] parallel libraries, 
is performed on an IBM SP2 parallel computer at the Pennsylvania State University Center for 
Academic Computing [10], Occupying four CPUs, the actual computation time is 4.45 sec- 
onds/iteration, which gives 7 microseconds/iteration/number of grid points. The computation 
performs 500,000 iterations before a fully developed wake flow is generated. 


B. Initial Triggering of Vortex Shedding 

In the early stages of the computation, it is difficult and time-consuming to initiate vortex 
shedding in the inviscid flow without some excitation. To start the first vortex shedding, the ve- 
locity normal to the base of the cone is excited in the form 
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u Sr, &, 0 = ew. Sin 2 ^ j cos(26>) exp[- (In 2)(/,r) 2 ] sin(2#» , (0 < r < , 0 < 6> < 2;r) 


(15) 


where £ is the amplitude of the fluctuation, is the free-stream velocity, r is the radial distance, 
R is the base radius, $ is the azimuthal angle, and f s is the frequency of the excitation. In Eq. (15), 
the distribution of the excitation velocity has four lobes with alternating signs on the cone base, 
and the amplitude decreases rapidly to half its initial value in one time period. In the present 
computation, the initial amplitude is chosen as e = 0.05 and the frequency is/ s = 0.171 u m /D. 
This corresponds to the regular vortex shedding frequency or the Strouhal number observed by 
Calvert [1], Soon after the excitation vanishes, a non-axisymmetric flow field develops and the 
first vortex shedding occurs. 

C. Investigation of Flow Field 

The unsteady wake flow may be visualized by the instantaneous pressure and Mach num- 
ber contours shown in Figs. 4 and 5 at the final time step (non-dimensional time is t' = uJ/D = 
64.04). Figure 4 shows the longitudinal distribution in a side view, and Fig. 5 shows the circum- 
ferential distribution in a rear view. The instantaneous three-dimensional unsteady flow field is 
very different from the two-dimensional axisymmetric steady flow field shown in Fig. 2. The 
pressure contours show many irregular vortices generated from the edge of the cone base and 
the Mach number contours show strong flow fluctuations behind the cone base. The flow pat- 
tern appears to be random or chaotic unlike the regular Kirmdn vortex street [11] behind a two- 
dimensional bluff body. It can be seen that some vortices pair and merge mutually in the down- 
stream region. This results in a large-scale motion of the far wake. However, the vortices very 
near the outflow boundary are dissipated non-physically because of the lack of grid resolution. 

In order to investigate the unsteady motion of the wake structure, time-traced snapshots of 
vorticity magnitude iso-surfaces are presented in Fig. 6. Figure 6 shows some vortex rings near 


10 



the cone base that have non-axisymmetric distorted shapes, which seems to induce the chaotic 
flow fluctuations in the near wake and the spontaneous vortex shedding downstream. It is also 
shown that the vortex rings in the near wake region interact each other and change their shapes 
into the axial vortex tubes in some transitional region away from the cone base at a distance 
nearly equal to the wall diameter ( x/D « 1). Further downstream, the vortex tubes shed from 
the transitional region possess a velocity component in the circumferential direction. This re- 
sults in a spiral motion of the far wake structure. Moreover, the vortex tubes sometimes twist, 
bend and stretch in the far wake region. These phenomena occur quite randomly and it is al- 
most impossible to pick up any two snapshots with the same instantaneous flow field during 
the entire computation. It is difficult to discern any periodicity from these time-traced observa- 
tions of the flow field. In the next subsection, a frequency analysis of the flow signals is pre- 
sented to quantify the periodicity. 

D. Analysis of Flow Signals 

The unsteady pressure and axial velocity along the centerline behind the cone are obtained 
as a function of time in the fully developed wake stage. By averaging the signals in time, the ax- 
ial variations of the mean pressure and axial velocity may be obtained. They are compared with 
Calvert's experimental data [1] in Figs. 7 and 8. Figure 7 shows that the mean pressure coeffi- 
cient from the present computation agrees with the experimental data well in the near wake re- 
gion up to the location of its minimum value. It then recovers faster to the free-stream value 
with a smaller overshoot than the experimental data in the far wake region. Figure 8 shows that 
the mean velocity curve of the present computation is a little more positive than the experimen- 
tal data and recovers slightly faster to the free-stream value too. The agreement between the 
predicted pressure coefficient and mean axial velocity is much better than that achieved in the 
low-order, unstructured grid calculations of Long et al [2]. The mean stagnation point, where 


11 



m = 0 and C p - 0 , of the present computation is slightly closer to the cone base than that of the 

experiment. It is likely that the small discrepancy between the present results and the experi- 
ment data is caused by the difference between the inviscid flow and the real viscous flow. The 
Euler equations have no physical viscosity that dissipates out the flow kinetic energy as the 
wakes move downstream, and the inviscid flow passes over a relatively smaller body in its ef- 
fective size than the real viscous flow since there is no boundary layer displacement effect. 

Two variables are defined in Ref. 1 to characterize the fluctuating components of the flow. 
They are denoted by, 

pr* — 

0 xlOO and <p=~ — xlOO 

u - K 

where u = w m — u m is the fluctuation of the axial velocity about its mean value and u m = |w| is 
the magnitude of the axial velocity. Thus, 9 is a measure of the absolute level of the velocity 
fluctuations and 0 represents the local turbulence intensity. The axial variations of 6 and pare 
shown and compared with the experimental data in Figs. 9 and 10, respectively. Figure 9 shows 
that the present results agree well overall with the experimental data except for a slight over- 
shoot in the near wake region and a small underprediction in the far wake region. Figure 10 
shows the same kind of discrepancy between the present results and the experimental data in 
terms of the shortening of the axial development as is explained for the mean flow results in 
Figs. 7 and 8 above. The predicted overall amplitude of the fluctuations agrees very well with 
the measurements. It is much better than the predictions by Long et al. [2] particularly in the far 
wake region. However, the present grid is much finer in that region and the order of accuracy 
of the present numerical scheme is higher. A comparison of Figs. 8 and 10 shows that the region 
of the highest local turbulence intensities is in the vicinity of the stagnation point. This is consis- 
tent with Calvert s [1] observations. However, this high value is associated with the lower mean 
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velocity at this location rather than an increase in the absolute value of the velocity fluctuations. 
In fact, as seen in Fig. 9, there is a slight decrease in the absolute value of the fluctuations in this 
region. 

The point of the mean pressure minimum in Fig. 7 corresponds to the transitional region 
(x/D * 1) where the axial vortex tubes start to form as described in the previous subsection. It 
coincides with the point of the mean velocity minimum where the highest speed of reversed 
flow on the centerline occurs as shown in Fig. 8. This point has significance in estimating the 
starting position of the periodic vortex shedding. Calvert [1] remarked that a periodic wake mo- 
tion first appears in the region of the mean pressure minimum and the periodicity presumably 
arises from the instability of the free shear layer in an adverse pressure gradient. He also men- 
tioned that the periodicity is most prominent in the region of the mean pressure maximum. In 
order to quantify the periodicity, the frequency spectrum of the axial velocity is obtained at the 
point of the mean pressure maximum (x/D = 2.0). This is shown in Fig. 11. The highest peak in 
Fig. 11 shows that there is a strong periodic wake motion at a non-dimensional frequency of 
fD/u„ = 0.171. This is exactly the Strouhal number found by Calvert [1]. The existence of this 
peak is taken as further evidence that the present computation succeeds in simulating the vor- 
tex shedding process accurately. 


V. Conclusions 

The subsonic inviscid wake flow past a cone has been simulated using high-order finite dif- 
ference schemes on structured grids. The present computation, performed in an axisymmetric 
structured grid system is achieved by changing the form of the flux vectors in the Euler equa- 
tions to remove the centerline singularity in the generalized coordinates. This approach makes 
it possible to investigate the complex wake flow field and to obtain the accurate values of the 
flow properties. It is shown that the vortex rings in the near wake change their shapes into axial 
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vortex tubes in the far wake through a transitional region where periodic vortex shedding be- 
gins. A spiral motion of the wake is found. The mean flow pattern agrees well with the experi- 
mental data when account is taken of the likely differences between the inviscid and the viscous 
cases. It is confirmed that the point of the mean pressure minimum is the starting position of the 
periodic vortex shedding. A spectral analysis of this periodic phenomenon shows that the 
Strouhal number is 0.171. This is in exact agreement with the experimental observation. On the 
basis of the very good agreement between the predictions and experiment, it is expected that 
the present methodology may be used for further analysis of wake-dominated flows past an ax- 
isymmetric blunt-based body. 
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Table Caption 

Table 1. List of optimized coefficients for compact finite difference schemes 

Figure Captions 

Fig. 1. Diagram of computational domain 

Fig. 2. Grid mesh system: (a) entire view, and (b) zoomed view 

Fig. 3. Zoomed side view of axisymmetric steady-state flow field for the initial condition: (a) 
pressure (p/p„), and (b) Mach number contours 

Fig. 4. Zoomed side view of instantaneous flow field at the final time step: (a) pressure (p /p_), 
and (b) Mach number contours 

Fig. 5. Zoomed rear view of instantaneous flow field at the final time step ( x/D = 1.0 from the 
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cone base): (a) pressure (p/p~), and (b) Mach number contours 
Fig. 6. Time-traced snapshots of vorticity magnitude iso-surfaces (17 levels of |<3| D/u m from 
0.0 to 0.002) in the interval of time -t„ = 0.642 D/u„ between two successive pictures 
Fig. 7. Axial variation of mean pressure along the centerline behind the cone 
Fig. 8. Axial variation of mean axial velocity along the centerline behind the cone 
Fig. 9. Axial variation of velocity fluctuations in terms of free-stream velocity 
Fig. 10. Axial variation of velocity fluctuations in terms of the local axial velocity 
Fig. 11. Spectrum of axial velocity fluctuations at the point of mean pressure maximum 
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Equation (8) 

Equation (9) 

Equation (10) 

Equation (11) 

fl] = 0.6511278808920836 

flo,i = -3.061503488555582 

*i,o = -0.5401943305881343 

fl 2 .o = -0.1327404414078232 

fl 2 = 0.2487500014377899 

ooj = 5.917946021057852 

flu - 0.8952361063034303 

02.2 = -0.6819452549637237 

«3 = 0.006144796612699781 

= 0.4176795271056629 

o u = 0.2553815577627246 

fl2.3 = 0.7109139355526556 

a =0.5775233202590945 

06.i = 5.870156099940824 

fli,4 = 0.007549029394582539 

on = 0.2459462758541114 

0 = 0.08953895334666784 

fa = 3.157271034936285 

ai,o = 0.1663921564068434 
oti.2 = 0.7162501763222718 
A.3 = 0.08619830787164529 

02,5 = 0.003965415751510620 
0i,o = 0.03447751898726934 
Oi,, = 0.4406854601950040 
cti3 = 0.6055509079866320 
0i A = 0.08141498512587530 


Table 1 




J. W. Kim and P. J. Morris 



Figure 1 




J. W. Kim and P. J. Morris 



Figure 2-(a) 
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Figure 2-(b) 
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Figure 3-(a) 
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Figure 4-(a) 
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Figure 4-(b) 
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Figure 5- (a) 
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Figure 5-(b) 
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Abstract 

The preliminary implementation and testing of the quadrature free Discontinuous 
Galerkin Finite Element Method in one and two dimensions are discussed. With an 
objective of understanding the method, algorithms to solve the linear scalar advection 
equation with periodic boundary conditions using the method have been written and 
executed successfully. A qualitative comparison of Lagrange polynomials and simple 
monomials as choices for the basis function set is made. Results from the algorithms 
using both Lagrange Polynomials and simple monomials as basis functions are also 
presented. 

1 Introduction to the Discontinuous Galerkin 
Method 

Convection dominated problems arise in applications as diverse as aeroacoustics, 
gas dynamics, meteorology, oceanography, turbulent flows, viscoelastic flows, magneto 
hydrodynamics and electromagnetism, among many others. Devising robust, accurate, 
and efficient numerical methods for solving these problems is not a trivial task for two 
reasons. The first is that the exact solution of nonlinear problems develop discontinu- 
ities after finite time, and the second is that these numerical methods might display 
complicated and erroneous solutions at these discontinuities. Thus, while developing 
these numerical methods, care must be taken to guarantee that the discontinuities of 
the approximate solutions are physically relevant, and that the appearance of discon- 
tinuities in the approximate solution does not induce spurious oscillations [1]. 

The Discontinuous Galerkin Finite Element Method (DG method) is one of those 
methods that are being developed to successfully address the issues mentioned above. 
The original DG method was developed in 1973, but it is only recently that enhance- 
ment and evolution of the method are taking place. The method has been proven to 
be well suited for high order accurate large time scale simulations. An important dis- 
tinction between the DG method and the conventional finite element methods is that 
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the resulting formulation is local to the element and the solution is not reconstructed 
by looking at the neighboring elements. Thus each element in the DG method can 
be thought of as an independent entity, merely requiring the boundary data from the 
surrounding elements. Since the DG method incorporates numerical fluxes and dis- 
continuous elements, it can be considered as a generalization of finite volume methods. 
Owing to its unique finite element nature, the DG method has numerous advantages 
over classical finite volume and finite element methods such as [1, 2]: 

• The DG method can be used to obtain uniformly high order accurate solutions. 

• The method is highly parallels able since the elements are discontinuous and the 
mass matrix is block diagonal and readily invertible. 

• It is suited for complex geometries since it can be used on unstructured grids. 
The method has also been shown to be immune to mesh discontinuities. It also 
requires simple treatment of the boundary conditions. 

• It can easily handle adaptive strategies, since refining or coarsening the grid can be 
achieved without considering the continuity restrictions typical to the conforming 
elements. 

• It has several useful mathematical properties with regards to stability and con- 
vergence. 

• The method is compact, since each element is independent. This compactness 
allows for a structured and simplified coding for the method. 

• The method allows for heterogeneity in the elements, that is, the order of accuracy, 
shape, and even the choice of governing equations can vary from element to 
element. 

Although the DG method is less susceptible to the problems that are common in 
finite difference schemes, it is not without a few weaknesses [2]. The method ha s been 
recognized as expensive, in terms of both computational operation count and storage 
requirements. Although theoretically the method can be applied to an element of any 
shape, the requirement of numerical quadrature for the integrals in the formulation 
has restricted the applications to hexahedral and quadrilateral elements. The recently 
developed [3] quadrature free DG method tries to resolve these problems by using 
polynomial basis functions, the product of which can be integrated without numerical 
quadrature. 

The rest of this report deals with the quadrature free DG method applied to 
the scalar advection equation, with periodic boundary conditions. Results for one- 
dimensional linear and a nonlinear advection equation are presented, followed by solu- 
tions of the two-dimensional linear advection equation on a square domain. 

2 Formulation of the DG Method 

Let the solution in an arbitrary domain be governed by a conservation equation of 
the form 


u t + V • F = 0 


( 1 ) 
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Let the domain be divided into smaller elements z) that span the domain. 

Let the solution in each element be approximated using an expansion of a basis set 
given by 


JV+l 

= bj u j ( 2 ) 

i= i 

B =>> {6 fc , 1 < k < N(p , d) + 1} (3) 

The number of basis functions N + 1 depends on the order of expansion denoted by 
p, and the number of dimensions, d. Application of the traditional Galerkin method 
to equation (1) using (2) and (3) gives, 

J bk (y>t + V • F^jdQ = 0 (4) 

for fc = 1 . . .N H- 1. 

Using (2) and integrating by parts, equation (4) can be recast as, 

j n h E b i (^J dn ~ j a V& * Fdn + J' b k F R ds = 0 (5) 

for k = 1 . . . N + 1. 

where ds is the normal surface area element (or the line element in 2D), and Fr 
is the Riemann flux vector through the surface element, which will be approximated 
using the flux values of the element and the neighboring element. 

It is convenient to carry out the integrations in the formulation (5) in transformed 
coordinates. In a conventional finite element method using quadrature, the integrals 
are calculated on a mapped element (denoted by A), but the final equation is evaluated 
and assembled in the real coordinates of f2. Since there is no assembling of the elements 
involved in the DG method, by mapping even the variables uj into the transformed 
coordinates, one can get a compact form of the equation (5). If the transformation from 
D to A is linear, the storage requirements will be minimized. This will be explained in 
detail below. 

Let the coordinates in the transformed space be £, rj and £. Let the unknown 
variables in the transformed element A (£,f?, C) be denoted as Vj. Then, 

y A b k E b 3 J l J l dA - J A J-'Vh FdA + J' b k F R ■ \J s \ds = 0 (6) 

for fc = 1 . . . N + 1, where J = » and J s is the Jacobian of the surface (or line) 

coordinate transformation. 

In particular, considering a linear advection equation in 2D, the flux is given by 

_ A ^ ^ ^ N+l N + 1 

Fn = ia x u + jbyU = F& = ia x v + jb y v = ia x bjVj + jb y ^ bjVj (7) 

j — i j=i 

and (6) can be written as, 
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where 

Mij = [ bjbklJldA ( 9 ) 

J A 

and 

K {j = / J^V&^IJIA (10) 

J A 

for k = 1 . . . N + 1, j = 1 . . . iV + 1. 

Due to the presence of the term J^ 1 , the matrix is made up of a sum of several 
matrices. The Riemann flux is approximated as a Lax-Friedrich’s flux of the form 

F R ds= i j(F f + F r )-n- a(V T - V t ) Jds (11) 

where Ft is the flux of the element to the left of the edge, and vecF s is that of the 
element to the right of the edge, n is the unit normal to the edge from left to right, 
and a is a smooth positive function. 


3 Computational Aspects of the DG Formula- 
tion 

The number of operations and the amount of storage required for the calculation 
of the integrals in the above formulation are determined by:"' 

1. the order of the coordinate transformation and 

2. the choice of basis functions. 

If the transformation is linear, then J is a matrix made up of constants and |J| 
is independent of £ or 77 . This means that a single mass matrix and component 
stiffness matrices K%j can be used for all the elements, thus minimizing the required 
storage. This is accomplished by using arbitrary triangles with straight edges for Q 
and an equilateral triangle for A, with the origin at its centroid. This is displayed in 
Fig. 1. 

{i -» } = ^ ^ } • ^ = 

where 

a.l = X3 — £2, a 2 = l/sqrt 3 ( 2 x\ - X2 — £3) 

W = V 3 ~ 2/2, b 2 = l/sqrt 3 ( 2 yi - y 2 - 2/3) 

where x \ ... £3 and y\ ... 1/3 are as shown in the figure 1, and xo and yo are the 
coordinates of the centroid of the element Q. The transformation matrix Jt will be 


al a 2 
bl 62 
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Figure 1: Sketch of the transformation from 


physical to canonical coordinates 


different from, but related, to the transformation Jacobian J which appears in Eqn. 

(8). ( J T = J*). 

This linear transformation also greatly reduces the computational time requirement. 
Since a single mass matrix is involved for all elements, it can be inverted prior to the 
main computation. Hence time marching for each element reduces from solving a 
system of equations to multiplying and adding matrices and vectors. 

3.1 Choice of the basis functions 

The choice of the expansion basis set greatly influences the computation time. Re- 
searchers [4] have worked on various polynomial basis function sets varying from simple 
monomials to orthogonal Legendre functions. A comparison of only the conventional 
Lagrange polynomials and simple monomials is presented here. 

For a problem involving nonlinear flux , the integrals in M {j and K tJ will involve 
products of more than two basis functions in their integrand. Choosing the conven- 
tional Lagrange polynomials for the basis set will necessitate excessive computation 
using numerical quadrature at each time step. Even if the integrals are computed an- 
alytically prior to the main computation, a huge storage will be required. The newly 
developed quadrature free approach aims at solving this problem by choosing simple 
monomial expansions for the basis set, thus making the evaluation of integrals simpler, 
and the matrices M and K sparser. The monomial expansions will be of the form 
fV* (For e.g. a second order basis set can be B = {1, 77, ^77, 77 s }. However, there 

are two drawbacks of the quadrature free approach. The first is that the condition 
number of the mass matrix is very high when monomials are used as basis functions. 
The second is that the evaluation of boundary integrals in Eqn. (8) will involve edge 
transformations. Hence, evaluating the boundary integrals in the quadrature free ap- 
proach will require relatively complex coding and more computational effort. But this 
increase in computation effort will be offset by the decrease owing to the simplified 
calculation of the mass and stiffness matrix integrals. 

On the other hand, if the flux is linear , using either Lagrange polynomials or the 
monomials, the integrals in M xj and K x j can be calculated analytically and stored 
before the main computation, without much usage of the memory. The mass matrix 
computed using Lagrange polynomials will be denser, but with a relatively lesser condi- 
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tion number compared to the mass matrix computed using monomials. The calculation 
of the boundary integral will be easier and faster using the Lagrange polynomials. 

4 Results From the Implementation of the DG 
Method 

Several programs have been written to implement and test the quadrature free DG 
method in one and two dimensions. The results from these programs for the linear 
scalar advection equation in one and two dimensions using both the Lagrange polyno- 
mials and monomials are presented. A nonlinear advection equation in one dimension 
is solved using the quadrature free approach. A periodic boundary condition is imple- 
mented for both one and two-dimensional calculations. Time marching is accomplished 
using Crank Nicolson’s 2 nd order scheme, a Fourth— order Runge— Kutta explicit scheme, 
and the Total Variation Bounded Runge-Kutta three (TVBRK3) stage method. 


4.1 One dimensional calculations 

The test equations are 


du du 

m +a di = 0 


( 12 ) 


and 


du 

m +a 


d{v?/2) 

dx 


= 0 


(13) 


on the domain 0 < x < 10 with periodic boundary conditions. 

The solution to the linear problem (12) is obtained with two initial conditions. The 
first is a half sine wave over a unit interval, and the second is a normal shock. While 
using the monomials, the initial condition was expanded as a Taylor series about the 
center of each element. The solution obtained using Lagrange polynomials is shown in 
figures 2a(l s£ order approximation polynomials with 2 nd order time marching) and 2b 
(2 nd order approximation polynomials with A th order time marching). 

The programs for the implementation of the DG method in one dimension are 
written in Matlab. The interval from x = (1, 10) is divided into 100 parts: that is, the 
value of Ax = 0.1. It is observed that the time step required for stability decreases as 
the order of spatial approximation increases. 

The calculations for the linear advection equation with monomials are shown in 
Fig. 3a, b and c. It was observed that the use of monomials lead to oscillations in 
the approximate solution where the exact solution had sharp discontinuities, as can be 
seen in Fig. 3b. Finally, the solution for the nonlinear Eqn. 13 is shown in Fig 3d. The 
initial condition was a sine wave, which transformed into an ‘N’— wave after sufficient 
time. 
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a 


b 


Figure 2: One dimensional calculations using Lagrange polynomials 



a b 



c d 

Figure 3: One dimensional calculations using monomials, (a) 3 rd order approximation mono- 
mials with 3 rd order time marching, (b) Absolute error in the calculation, (c) Calculations 
with an initial condition of a normal shock wave, (d) Solution for the nonlinear advection 
equation. 
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a b 

Figure 4: The initial condition, (a) The sine wave pulse with unit maximum value at the 
center of the domain, (b) a contour map of the initial condition. 


4.2 Two dimensional calculations 

4.2.1 The advection equation 

The test equation is 


du 

dt 


“I" CL 


du du 


= 0 


(14) 


u(0,x,2/) = [sin(7rx) sin(7n/)] 4 

on a periodic square domain -1 < x, y < 1 . The grid is made up of 800 similar 
right angle triangles. The program is written in Fortran 90. The basis set is made up 
of Lagrange polynomials of up to second order and monomials of up to 4 th order. The 
initial condition is shown in Fig. 4. 

Figure 5 shows the results for the case with the basis set made up of linear monomi- 
als. A similar result is obtained with all other choices of basis functions. The speeds of 
the wave in the x and y directions are taken to be the same. It is observed that the size 
of the time step required for a stable solution decreases as the order of approximation 
is increased. 

4.3 Calculations for the linearized Euler equations in two 
dimensions 

Consider the propagation of an acoustic pulse in a constant mean flow. The lin- 
earized Euler equations are given by, 
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( 15 ) 


dU dE dF_ 
dt + dx dy 


with 




✓ \ 
p 


M x p + p 0 u 


MyP+POV 

u 

>,E=< 

M x u + p/p 0 

> P = i 

MyU 

V 

M x v 


MyV + p/p 0 

, p , 


M x p + p 0 u t 


„ MyP + PQV d 


where p,u,v,p are the flow variables, and M x and M y axe the values of the mean 
flow in the x and y directions respectively, non-dimensionalized by the speed of sound 
based on the mean thermodynamic values, and po is the mean flow density. 

The finite element discretization of the above differential equation results in a sys- 
tem of equations that can be written as, 


-M 

0 

0 

0 ' 


\ K n 

K n 

K 13 

0 ’ 


Fl 1 

0 

0 

M 

0 

0 

M 

0 

0 

dU 
~dt + 

0 

0 

K 22 
0 

0 

K 33 

Ki 4 
K 34 

U= ■ 

F 2 1 

. 0 

0 

0 

M. 


. 0 

K 42 

K 43 

K-u - 


F* J 


(16) 


where M and Kij are the mass and stiffness matrices, and F\ ... F 4 are the right 
hand side vectors of the corresponding individual equations. 


4.3.1 Results 

The domain is divided with unstructured triangulation. Basis sets of degree 1 
and 2 (order 2 and 3 respectively) are used. A square domain is considered, with a 
uniform mean flow from the left to right, over a rigid bottom wall. An initial Gaussian 
distributed acoustic pulse is located near the lower boundary, as shown in Fig 6. 

Reflective boundary conditions are imposed on the lower boundary, and non-reflecting 
boundary conditions are imposed on the three other boundaries. The radiating bound- 
ary conditions are implemented using characteristic variables [5]. Figures 7a and b 
show the computed propagation of the acoustic disturbance with the mean flow. Some 
reflections are observed near the right boundary, and these are due to the characteristic 
method of boundary condition used, as mentioned in [5]. 


5 Conclusions and Future Work 

The linear advection equation with periodic boundary conditions has been solved 
using the DG method in one and two dimensions. A qualitative comparison of the 
Lagrange polynomials and simple monomials as choices for the basis set has been made. 
A basis set comprised of simple monomials is more suited for h-p refinement of the 
solution. The choice of monomials for nonlinear convective equations is thought to be 
advantageous compared to the Lagrange polynomials under the formulation described, 
in which, the nonlinear flux is evaluated by multiplying or dividing the approximating 
expressions for the basic variables. However, the formulation involving the nonlinear 
flux could also be solved using iterative techniques. It remains to be investigated 
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